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UNSTEADY NATURAL CONVECTION IN CYLINDRICAL ENCLOSURES 
WITH STRATIFIED MEDIUM 


By 
RUEY-JONG SHY 
August, 1985 
Chairman: C.K. Hsien 
Major Department: Mechanical Engineering 
A numerical study was given to investigate unsteady natural 
convection in enclosures with thermally stratified water. Two 
different designs of enclosure wall were analyzed, one placing 
insulation on the outside of the wall, the other placing the 
insulation over the inside of the wall. An uninsulated enclosure was 
also studied for use as control. For all the cases studied, water 
was stratified at the midsection of the enclosure; hot water occupied 
the top section, while cold water was located at the bottom 
section. The wall temperature was initially uniform at the hot-water 
temperature. 
The problem was solved as a conjugate heat transfer problem to 
determine the flow and temperature fields in water ana the 
temperature field in wall. in the numerical solution, the wall and 


‘the water were treated as two separate regions, and conditions were 


cy 


derived for matching temperature and heat flux at the wall/water 
interface. The governing equations were transformed into a 
vorticity-stream function formulation, and control-volume approach 
Was Used to derive the Fini te-difference equations. The wall and 
water energy equations and the vorticity transport equation were 
solved by the use of Alternating Direction Implicit method, while the 
stream function equation was solved by the use of Successive Over- 
Relaxation method, 

The numerical results indicate that exterior and interior 
insulation are about equal in preserving energy; however, the 
interior insulation performs far better in retaining stratification 
in water, an outcome that is highly desirable for solar/therma) 
energy-storage applications. The velocity and temperature fields in 
‘the boundary layer at the early stage of cooling were also studied, 
and physical mechanisms leading to temperature defect and flow 
reversal are presented. The numerical results were checked by 
‘vesting stability and convergence, and they were also compared with 


experimental data reported in the literature. 
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CHAPTER I 
INTRODUCTION 


Natural convection problems are usually divided into to 
categories: external problens, such as flow and heat transfer fron a 
heated rod or plate to a fluid at rest; and internal problems, such 
as flow and heat transfer between parallel plates or ina fluid= 
filled enclosure. Much of the earlier work on convection studies tas 
been concentrated on the external problems. The internal provlens 
have received attention nore recently and it is due to the fact that 
the convection in an enclosed space is complicated. 

For an external problem of natural convection from a heated rod 
or plate toa fluid at rest, the Prandt! boundary layer theory can be 
used to simplify its governing equations. The ambient condition is 
also unaffected by the heat transfer. For an internal problem of 
natural convection in a confined space, however, a core region is 
surrounded by the boundary-layer regions near the walls, and the 
coupling of these two regions leads to a great complexity in heat 
transfer analysis, 

The internal convection problems have extensive applications in 
geophysics and engineering. For example, geophysicists are 
interested in studying the Fluid motion induced by density gradient 
in a cavity. Solar energy specialists are interested in suppressing 


the natural-convection heat loss in solar collectors. With the 


recent miniaturization of electronic equipment, electrical engineers 
also become interested in the convection problems. Certainly, these 
convection proplems remain a major area of research to heat transfer 
engineers, and the recently developed numerical methods enable them 
‘to model complicated convection problems more realistically; as 4 
result, there is an outburst of convection research as reported in 
the literature (a review follows in the next chapter). 

It is the purpose of this dissertation to study numerically the 
natural convection of water in a stratified enclosure. The nunerical 
solution is applied to compare two alternatives in the design of 
solar energy storage tanks, one placing insulation on the outside 
surface of the tank (a common practice), and the other placing the 
insulation over the inside of the tank wall (a new idea), Since it 
is important to maintain stratification in the storage tank in order 
to raise the overall efficiency of a solar energy system, it will be 
shown that significant improvement in stratification can be realized 
by placing the insulation over the interior of the tank wall. It is 
hoped that the insight gained from this research will be useful in 


the future design of solar energy storage systems. 


CHAPTER II 
LITERATURE. REVIEW 


The internal convection problems have been extensively reviewed 
by Ustrach [1] in 1972 and by Catton [2] in 1978. Since the present 
dissertation will deal with a study of unsteady natural convection in 
enclosures with stratified medium, and the study is related to the 

analysis of storage of thermal eneray for solar energy applications, 
‘the convection review that is to follow will address these tio areas 


of research as docunented in the literature. 


2.1 Natural Convection in Enclosures 
For natural convection in rectangular and cylindrical enclosures 
bounded by vertical side walls (i.e., vertical enclosures) either 
heated or cooled, two types of problems usually arise and are 
reviewed here. The first type of problem deals with side walls 
‘imposed with different temperatures, commonly known as differentially 
heated conditions; see Figure 2.1(a). The second type of problem has 
these walls imposed with symmetrical conditions, either constant 
‘temperature or constant heat flux; for ease of reference later to 
this second problem, these will be called symmetrically heated 
Conditions; see an example shown in Figure 2.1(b). The top and 


bottom boundaries for both problems are usually insulated. 
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Figure 2.1. Natural convection in enclosures 


2.1.1 Differentially Heated Conditions 


In 1954, Batchelor [3] studied natural convection in gas inside 
a rectangular enclosure imposed with isothermal side walls and either 
Perfectly conducting or insulating top and bottom surfaces. From the 
analytical results obtained, he reasoned that several flow regimes 
occur inside the enclosure depending on the aspect ratio of the 
enclosure and the Rayleigh number. Eckert and Carlson [4] tested 
with enclosures made up of isothermal side walls at temperature 
differences ranging from 5* to 90°C. They were able to define 
Conduction, transition, and boundary-layer regimes as functions of 
‘the Grashof nunber and the aspect ratio, The enclosed core region, 
however, was not isothermal as predicted by gatchelor, but stratified 
‘in that the temperature was uniform horizontally but increased in a 
vertical direction, Elder [5] made a study of natural convection in 
water in a vertical slot with aspect ratios ranging from 1 to 60 and 
obtained basically the same results as Eckert and Carlson. An 
analytical study by Gi11 [6] provided further proof of Eckert and 
Carison's observations. 

The convection research carried out so far has established that 
‘the basic flows inside an enclosure consist of two boundary layers 
near the heated or cooled walls, and the fluid in the enclosed core 
region is temperature stratified; see Figure 2.1(a). Schinkel et al. 
(7] studied the effect of core stratification on natural convection 
and showed that the stratification in the core region had an effect 
of decreasing ouoyancy force and producing secondary flows in the 


fluid. ohn et al. [8] were able to show that the boundary-layer 


Flows near the side walls were actually responsible for the heat 
transfer across the enclosure. For the boundary condi tions 
iMustrated in Figure 2.1(a), the fluid moves up near the heated 
wall, drifts across the top of the enclosure, and then descends along 
the cooled wall at the other side. 

The construction of the enclosure wall also has an effect on the 
convective flow in the enclosure. Gdalevich et al. [9] studied 
numerically the heat transfer ina vertically bounded air space; only 
one of the side walls had a finite thickness and the thermal 
conductivity of this wall was treated as a variable, the outside 
surface of this wall being imposed with an isothermal condition. The 
other side wal] was also isothermal but at a different temperature. 
As expected, an increase of the thermal resistance in the thick wall 
Jeads to a drop in the intensity of the convective motion in the 
enclosure. 

In the event that the ‘top and bottom end walls of the enclosure 
are not adiabatic, then the constructions of these walls also have a 
bearing on the heat transfer. Meyer et al. [10] have shown that the 
convective coefficients for enclosures fitted with adiabatic end 
walls are substantially higher than those with conducting walls. A 
study by Elsherbiny et al. [11] gave qualitative results for wall 
conduction, radiation, and convection interactions; they also 
Proposed a list of parameters that must be specified when reporting 
data for enclosures bounded by conducting end walls. More recently, 
Kim and Viskanta (12) reported both experimental and numerical 


results on the effects of wall conductance on natural convection in a 


rectangular cavity. They tested three different heating and cooling 
arrangements in which the external wall of the cavity was heated from 
either the side, the top or the bottom, while at the same time, the 
heat was removed from the opposing wall of the cavity. They showed 
‘that the neat conduction in the connecting (unheated) walls nad an 
effect of either stabilizing or destabilizing the fluid motion in the 
cavity, 

Most of the convection problems studied in the literature deal 
with two-dimensional systems in which the enclosure is treated as 
infinitely long in the third dimension. In reality, however, the 
enclosure is of finite length, and the third dimension is expected to 
have an effect on the fluid motion. Morrison and Tran (13] used a 
laser Doppler anemone ter to measure velocity profiles and showed that 
three-dimensional flow was induced by the heat loss from the vertical 
end walls in the third dimension. Horizontal velocity component 
formal to the end walls is small as compared with the vertical 
Component when these walls are well insulated, However, when they 
are not, as those found in many heat transfer experiments, a strong 
three-dimensional motion is detected. Mallinson and Vahl Davis (14) 
used a numerical method to show that the end effect was mainly 
confined to the vicinity of the end walls for air or high Prandtl 
fluids. A two-dimensional model is considered adequate for 
Convection analysis as long as the Rayleigh and Prandtl numbers are 


large and the cavity is long in the third dimension. 


2.1.2 Synmetrically Heated Conditions 


When the enclosure is symmetrically heated or cooled from the 
side walls, top and bottom surfaces being insulated as shown in 
Figure 2.1(b), the core region is again surrounded by the boundary 
regions. Schwind and Vliet [15] used both the Schlieren and 
shadowgraph techniques to study the natural convection and 
Stratification of water ina vessel with heated side walls and a free 
‘top water surface. The fluid motion leading to the formation of 
stratified layers was studied in great detail, and a reverse shear 
flow was found near the boundary layers in the presence of a large 
adverse temperature gradient in the core region. Evans et al. [16] 
tested with liquids in a partially filled cylindrical container with 
uniform heat flux applied on the side wall. They developed an 
analytical model for studying the fluid flow in the container. To 
simplify their analysis, they divided the container into three 
regions: a thin boundary layer near the side wall, a mixing region 
at the top where the boundary layer discharges and mixes with the 
upper core fluid, and a main core region which gradually diminishes 
in height because of the mixing of the fluid from above. Numerical 
results computed with this model were found in close agreement with 
experiment. Evans’ work was further extended by Hess and Miller (17) 
to cover heat transfer in the container wall. A thermocline-type 
cylindrical enclosure with its side wall imposed with an isothermal 
condition was solved with an integral method, and the results were 


found in good agreement with the finite-difference solution and the 


experimental data which Hess and Miller obtained in an earlier study 
ts). 

Transient natural convection in a cavity has also been an area 
of interest. darakat and Clark [19] made a numerical study of 
natural convection in a partially filled liquid container whose walls 
were subjected to arbitrary time and space-dependent tempera ture 
variations. An explicit finite-difference technique was developed 
for the solution of this problem. A method of stability analysis 
Suitable for solving partial differential equations with variable 
coefficients was also included in their paper. An interesting study 
was made by Sakurai and Matsuda [20] that involved Boussinesq fluid 
ina circular cylinder. The top of this cylinder was ata higher 
temperature than the bottom so the temperature distribution in the 
fluid was thermally stable. The side wall was made up of a good 
‘thermal conductor with small time constant, Ata certain moment, the 
top and bottom temperatures of this cylinder were changed suddenly by 
an equal amount but in opposite sign, and the convection in the 
cylinder was studied as it responded to the temperature change in the 
side wall. They were able to show analytically that a meridional 
Circulation in the core region was induced by the boundary-layer flow 
near the side wall. This circulation redistributed the fluid to the 
extent that a new state of stratification took place in the core 
region. Sakurai and Matsuda's work was further extended by Jischke 
and Doty [21] who studied the transient flow in an arbitrarily-shaped 
container, completely filled with fluid. The temperature of the 


container wall was impulsively changed and it was found that the 


10 


invisid interior region heated up by a convective process that was 
riven by suction induced by the boundary layer. The inviscid, 
adiabatic interior, however, responded to a horizontal average 
temperature perturbation of the wall because the boundary layer 
smeared out the circumferential temperature variation. The interior 
thus, in effect, responded to an isothermal boundary in each 
horizontal plane. Four examples were provided in their paper for 
substantiation of their theory. 

Hess and Miller [18] used a laser Ooppler velocimeter to measure 
velocity field inside a thermocline-type cylindrical enclosure. Both 
the axial and radial components of the velocity were measured and 
found to have a magnitude ranging from 0.01 to 0.45 cn/s. Good 


agreement was also found with a numerical solution. 


2.2 Energy Storage in a Thermocline or a Stratified Medium 


Storage of solar energy in a stratified medium is always desir- 
able in air-heating systems and domestic hot water systems [22]. In 
practice, the heated medium from a solar collector enters the top of 
a stratified storage tank, and the medium leaving the botton of the 
tank enters the collector at a lower temperature. This arrangement 
permits the collector to be operated at a constant, high efficiency, 
which is important to the performance of the total system. 

Numerous models have been developed to simulate the fluid flow 
in a stratified tank. Close [23] developed an analytical model in 
which the tank was divided into several well-mixed sections. By 
postulating that the incoming hot fluid would seek a section for home 


that has the same density level as the incoming fluid, he was able to 


a 


develop partial differential equations for energy balance in each 
section, A three-section modal was found adequate for simulating the 
performance of a water heating system [24]. Close's analytical model 
is considered useful when the fluid velocity entering the storage 
tank is small, When this velocity is large, however, turbulent 
mixing effect must be accounted for. Hu and Han [25] considered 
viscous entrainment in their mathematical model. A set of partial 
differential equations was developed by them and solved with an 
implicit finite-difference method. 

Cabel1i [26] developed a two-dimensional model for studying the 
fluid flow and temperature distribution ina storage tank. Goth 
charging and discharging situations were modelled numerically and it 
was found that a one-dimensional analysis, such as that developed by 
Close, should provide reasonable accuracy in predicting the 
temperature profiles in a tank if the hot medium enters at the top of 


‘the tank. Sha and Lin (27] were able to develop a three~< 


imensional 
model for analyzing the flow stratification ina thermocline storage 
tank. They developed a computer program (CONMIX-SA) for analysis and 
‘the model was shown to give good prediction of the fluid Flow and 
‘temperature distribution in the tank during a demand mode. 

The benefit of thermal stratification has also been documented 
in the literature. Sharp and Loehrke [28] studied numerically the 
potential benefits of thermal stratification in a solar space heating 
application. They compared the performance of a system with 
stratified storage to that of a similar system with well-mixed 


storage and found that the benefit due to stratification, as measured 


R 


by the increase in the fraction of the total heating load supplied by 
solar energy, ranged from 5 to 15%. The benefit of thermal 
Stratification was also verified by Veltkamp [29] using both 
experiment and numerical modelling. For low temperature applications 
under a low flow-rate condition, a stratified storage was able to 
enhance the total system efficiency by 10% or more. They also noted 
that, in order to numerically model the storage tank with accuracy, 
itwas necessary to divide the tank into 20 sections at least, a 
proposal in sharp contradiction to the division recommendation made 
by Gutierrez et al. [24] in an earlier study. Cole and Bellinger 
[30] proposed a new strategy for controlling stratification that 
enabled lowering of the average collector operating temperature by as 
uch as 20°C, This reduction of collector temperature also resulted 
in an improvement in the total system performance up to 20%. 

The decay of thermal stratification has been studied in the 
literature, Viskanta and Hale [31] studied unsteady natural 
Convection ina thermally stratified tank driven by heat losses 
through tank walls. A two-dimensional model was used in their 
analysis and the results indicated that the natural convection was 
not confined to the vicinity of the walls but spread out to fill up 
the entire tank. To reduce the degradation of the thermocline and 
unwanted heat transfer between the hot and cold regions, the tank 
should be well insulated. Jaluria and Gupta [32] developed a one- 
dimensional model to analyze the decay of thermal stratification in a 
water body under a static condition. For the plexiglass container 


used in their analysis, they noted that the cooling process in the 


a 


water depended strongly on its initial condition. A stratified water 
body was able to curb mixing, and the energy transfer inside the 
stratified region was mainly a result of thermal diffusion. Gross 
[33] made a similar study of decay of stratification in a thermocline 
system, He determined the relative importance of each of the five 
heat loss mechanisms present in the system that include (1) heat loss 
to ambient, (2) convective motion induced by heat transfer from the 
hot layer to the cold layer via the highly conductive tank wall, (3) 
mixing during charging and discharging, (4) thermosiphoning loops in 
the inlet and outlet manifolds, and (5) thermal diffusion from the 
hot layer to the cold layer. Strategies to cut down on these losses 
were also proposed in his paper. 

Lavan and Thompson [34] studied experimentally the thermal 
Stratification in a hot-water storage tank under a high extraction~ 
rate condition. Parameters studied included the aspect ratio on the 
tank, the water inlet-outlet temperature difference, the inlet and 
exit port configuration, and the water mass flow rate. Sherman et 
a1. [35] made an experimental study of the effect of wall conductance 
fon the temperature relaxation time of a stratified tank under a 
Static condition. From the five sets of data they collected, they 
were able to show that conduction along the tank wall could cause 
decay of thermal stratification to a significant extent. Miller [36] 
showed that the natural convection generated at the fluid/wall 
interface of a stratified tank increased for larger Rayleigh numbers, 
but decreased for larger aspect ratios of the tank. This natural 


Convection also increased proportionally to the square root of the 


uu 


thickness of the tank wall. Chaney et al. [37] suggested the idea of 
insulating the interior surface of the storage-tank wall in order to 
cut down on the heat leak; however, they provided no verification of 


the usefulness of this idea. 


CHAPTER II 
MATHEMATICAL FORMULATION 


The system under investigation is shown schematically in Figure 
3.1, A cylindrical tank filled with water, with top and bottom 
surfaces insulated, is chosen for analysis. The water tenpera ture 
changes with time due to heat loss through the tank wall. In order 
to compare the performances of the tank with different designs of the 
wall, the insulation on the wall is placed either on the outside or 
on the inside of the wall, An uninsulated tank is also tested for 
control. The entire system analyzed 4s divided into three regions: 
the water occupies the space R;>r20, the wall/insulation is in the 
region RVR}, and the air is in roRy. 

The water is assumed initially stratified as shown by the 
isotherms plotted in the figure. The water bodies in the top and 
bottom sections of the tank are isothermal but at different 
(dimensiontess) temperatures. In the middle section, the water 
temperature is linearly stratified as shown by those horizontal Vines 
drawn, The initial temperature of the tank wall is assumed unity 
throughout, a condition necessary to provide perturbation in water to 
initiate convection. 

For the problem under investigation, the water is considered 
Viscous and having constant properties (except density); the natural 


convection is laminar. In the tank wall, the thermal diffusivity is 
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Figure 3.1 Coordinate system and initial temperature distribution 


cu 


treated as constant. The temperature distributions in the water and 
wall are to be investigated, and they are found by solving governing 


equations in these two regions. 


3.1 Governing Equations and Conditions for 


Natural Convection in Kater 


The conservation of mass, momentum, and energy in water may be 
expressed in three equations as follows [38]: 


Continuity Equation: 


Tet ered x0 ay 


Momentum Equa tion: 


ope F- me wy ek ony (3.2) 
Energy Equa tion: 


OT eye op age 
oC, DET WRIT + oT TP + ue + q (3.3) 

With the assumptions made earlier, the pressure, dissipation, and 

source terms in the energy equation my be neglected. In the 

momen tum equa tion, 


F- w= (p - 9,15 - 94 (3.4) 
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Boussinesq approximations are used to relate the density difference 


to temperature difference as 
p> 9, = BT, = T) (3.8) 


The density differential term in the continuity equation may be 


dropped, thus giving a new set of simplified equations: 


(3.6) 
OF Los 
odes PRAT, - T) = wy + wy (3.7) 
or 
aC, Be = kT (3.8) 


AS expected, the momentum and energy equations are coupled through 
‘the temperature terms in these equations. 

Equations (3.6) through (3.8) are specialized to the case of a 
‘two-dimensional cylindrical system of interest in this study. A set 
of four equations is written as follows: 


Continuity Equation: 


%y, . 
yd Za 
wrth tgito (3.9) 


Momentum Equation in r Direction: 


a 
v v, ay, ap, av 


rn esa Ma, te 
a ee ee oe) 
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ey, 
+ J (3.10) 
ae 
3. 2", wv. a. 2». 
2 2 Pe al es ae 
a ee eae tea 
ay, 
+ Fl + gBtT - 1) (3a) 
a 
Energy Equa tion: 
2, 
Tey Dey Me et, 1ar at 
ae* % * iz bee tr ae. (3,12) 


Eliminating pressure terms in the two momentum equations given 
above and relating the stream function, %, and the velocity 


components 


Lay 
= 


(3.13) 


ay 
rar 


and further introducing the vorticity relation 


(3.14) 
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give the vorticity transport equation 


2 
H+ vo) + ve) =e Te ed +P) (3.15) 
2 


The energy equation may also be recast as 


aT ,1a a waetad 
Fete H DY + Hn = ah et (3.16) 
It 4s comon practice in the fluid flow and heat transfer 
analysis to introduce dimensionless groups 
ig 
wed 
Per «8 
fo kt vf oth 
* F 
(3.17) 
agt 
ei, 
Car " 
oR 
pect +28 
os we 
% a) 


so that the number of variables to be considered in the analysis my 
be reduced. For the sake of convenience in expressions, the 
asterisks in those dimensionless quantities are dropped. Then, the 
vorticity transport equation, the stream function equation, the 


velocity equations, and the energy equation may be expressed as 
follows: 


a 
Vor ticity Transport Equation: 


2 
va) + Z4uw) = rca eid 2 ry) ¢ Ley 
% alte wpe) +e 


aT, 
Pe at (3.18) 


Stream Function Equation: 


1a aw. 
ee! +A wets ae w (3,19) 
Velocity Equations: 
ees 
wed. (3.20) 
rar; 
os: 
woes (3.21) 


£13 


(3.22) 


Te a ete 
rete + UT)» ACRE ae Seg ghd ed Leng 


The initial and boundary conditions for water may be summarized 


as 


Tylgs2.0) = Ty orgsz) 


Ey 


a, 
zel, geo, ogre. (3.23) 
at 
f 
ret0, 920, weo, wo 


reed, veo, Bao, rez) = tle 
Notice that there are five unknowns (us ys Uy vy and Tp) in the 
governing equations (3.18) through (3.22), and there are five 


equations here to solve them, 


3.2, Governing Equations and Conditions for 
aenetyen Conduction ta aT 
The energy equation in the tank wall is the heat conduction 
equa tion 


2, 
He Lig Maly , a20gt) 


He parr rs (3.24) 


in which the diffusivity, a, has been placed after the partial 
differential operators in order to allow its use in the analysis of 
heat conduction in walls with or without insulation. 


The following dimensionless groups are introduced for the wall 


rR 


(3.28) 
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and once again, the asterisks are dropped and equation (3.24) is 


recast as 
at, star) aat) — a%(aty 
ig ae a Pure s 
ee ee 
re 2 
where 
gmc 
5 
and 
AR, 
Bs 
ae 


(3.26) 


Initial and boundary conditions for the wall are given by the 


relations 


Tylrg.2,0) = Telrgyz) 


aT, 

z=0, x0, 
aT, 

zl, so, 


=O, TelOyz,0) = Ty(z.t) 


(3.27) 
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Notice that the dimensionless radius in the wall has been redefined 
with reference to the wall itself; such a change permits the 
conditions of the wall to be expressed mare conpactly as those given 
above. 

For the problem given, only the anbient air temperature, T,, is 
fixed; the water temperature is constantly changing with tine due to 
hneat loss through the tank wall. In fact, the wall temperature also 
changes, but the heat fluxes and temperatures at the wall/water 
interface renain continuous at any instant of time, This situation 
Presents a conjugate heat transfer problem in that conduction in the 
wall and convection in the water must be considered simul taneously, 
with proper matching of conditions at the interface. This problem 
can only be solved numerically as will be shown in the next two 


chapters. 


CHAPTER IV 
METHODS OF SOLUTION 


Several problems are encountered and mist be resolved before the 
governing equations given in Chapter III can be solved numerically. 
Also, these governing (differential) equations must be discretized 
and appropriate numerical methods must be chosen to solve them. A117 


‘these must be done at the early stage of the solution of the problem. 


4.1 Special Treatment in the Solution 


Special treatment is needed to resolve several problems 
described in this section. 


4.1.1 Computation of Heat Loss from Wall to the Ambient Air 


It must be noted that, in the formulation of the governing 
equations given in the preceding chapter, only those in the storage 
medium (water) and the tank wall have been given. The external 
convection from the wall to the ambient air has not been considered 
a5 a separate region but accounted for in the boundary conditions for 
the wall, Certainly, it is possible to consider the external 
convection as a separate problem and solve it together with the wall 
and water problems; however, such a decision would require matching 
of extra conditions at the wall/air interface in addition to the 
solution of the heat transfer in the air, an exceedingly laborious 


task even for a computer to handle. The present treatment of the 


a 
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external convection by means of use of proper boundary conditions for 
the wall is therefore preferred, 

In the search of empirical equations for estimating wall heat 
loss, attention must be given to the fact that the wall is 
Nonisothermal. In reality, the outside wall temperature is the 
result of the wall heat flux, which is, in turn, related to the 
temperature of the wall itself. It is thus desirable to find 
empirical relations that permit a simultaneous determination of the 
‘temperature and heat flux in the wall. 

The external convection from a nonisothermal, vertical wall has 
been a subject of numerous studies (39-43). Of these references 
cited, Sparrow and Gregg [39] and karvinen's analyses [40] are chosen 
for comparison. Sparrow and Gregg used a similarity analysis to 
evaluate the Nusselt number for natural convection from a vertical 
wall on which the temperature distribution followed a power-law 
relationship 


Ty - Teo He a 


Sparrow and Gregg presented their results graphically; a fit of their 


curves using a cubic polynomial yields the relation 


Nu 
—_ = 0.02568 m? - 0,135 n® + 0.3523 m + 0.5018 (4.2) 
Ds 


Gr 
GR 


Notice that, although the temperature of the wall varies with N in 


2 


equation (4.1), the Nusselt-nunber equation correlated above is 
independent of it. 

Karvinen used the approximate solution methods of Raithby and 
Hollands [44] and developed a relation between the wall heat flux and 
wall temperature distribution 


iS i Ak 
BF ar cy° ar az (4.3) 
° 


az) = ck¢ 


where Prandt] nunber, Pr, is greater than 0.6; 47, stands for a 


temperature difference, T, - Te. The constant C is related to the 


Prandtl number as follows: 


Cs 8.503 (4.4) 


eT 
ta + 882° 


By modelling the temperature variation on the wall as a series of 


linear functions, Karvinen was able to integrate equation (4.3) to 


sive 
81 + = e $ 3 -4 
Pr 3 
oh a oe eet (4.5) 
where 
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Equation (4.5) may be used together with Newton's cooling law to 
evaluate the local convective coefficient. 
For a power-law temperature variation given as equation (4.1), 


equation (4.3) may be integrated to give 


a 34 P 
Ey 0.4986 eta (4.8) 
or, * 
me) 


This equation may be used to test the accuracy of the approximate 
solution by Karvinen, Numerical results are obtained for both 
solution methods in which the Prandtl nunber is taken to be 0.7 
(air). AS 45 shown in Table 4.1, the error in the Nu,/(Gr,/4) YA 
values computed with the approximate solution increases with m; for 
all the m values tested, the approximate solution always gives a 
lower estimation of the WusseTt nunder. 

In selecting a prediction method from the to approaches 
compared above, consideration was given to the fact that (1) it is 
difficult to fit the wall temperature using the temperature relation 
given by equation (4.1), and (2) even for situations where the wall 
temperature can be fitted by this equation, the exponent m for z is 
usually smal] (< 0.2), For these reasons, Karvinen's equation (4.5) 
is preferred. A point should be noted, however. The fact that his 
equation was derived on the basis of heat transfer froma flat wall 


does not affect its use in the present analysis because the storage 
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Table 4.1 Comparison of exact solution and approximate solution for 
heat loss calculation 


YN 
Nu, /(6r/4) 
m | exact (eq. (4.2)) | approximate (eq. (4.6))] error & 
o 0.5018 0.4986 0.64 
0.1 0.5357 0.5182 3.27 
0.2 0.5671 0.5358 5.52 
0.3 0.5960 0.5518 7a? 
4 0.6228 0.5665 9.04 
0.5 0.6474 0.5802 10.38 
a oa 0.5929 11.52 
07 0.6911 0.6089 12.47 
0.8 0.7104 0.6162 13.26 
0.9 0.7282 0.6269 13.31 
1.0 0.7448 0.0872 14.45 


* isothermal condition 
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‘tank under investigation has a small curvature, according to a 


numerical study made by Cebeci [45]. 


4.1.2 Hatehing of Conditions at vall/Water and al/tnsulacion 
‘Tnterfaces 


For the conjugate heat transfer problem solved in this work, 
both the heat fluxes and the temperatures must be continuous at all 
interfaces. As has been mentioned before, the heat flux is related 
to the temperature, neither of these two quantities, therefore, being 
known a priori. 


The interface conditions may be matched by different methods as 


Feported in the literature; see references in the literature 
(36,46,47, for example]. The method used in this work to match 
wall/water interface conditions is derived on the basis of the use of 


continuity of temperature 
T,(0.2,t) = Te(1,2,t) = T, (zt) (4.7) 


and heat flux 


aT, aT, 


(4.8) 


By writing equation (4.8) ina finite-difference form and 
incorporating it with (4.7), it is possible to write the interface 


‘temperature a5 
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(4.9) 


where Tp.) and Ts; represent temperatures one grid point avay from 
the interface in the fluid and the solid region, respectively; Sr 
and rg refer to unit lengths in these regions. In the numerical 
work to be presented later, this interface temperature, Tj, 1s used 
4S a conmon boundary condition for both wall and water regions, and 
the wall and water energy equations are solved simultaneously. Heat 
fluxes across the interface are then computed, and if they fail to 
agree, the newly computed temperatures are used to recalculate 
interface temperature using equation (4.9). This procedure is 
repeated until all heat fluxes match at the interface. 

For the interface between the wall and the insulation, a 
different method is used. By considering the wall and the insulation 
as a single region in the numerical computation, it is possible to 
use Patankar's contro! volune approach [46] to lay out the grid 
points as shown in Figure 4.1. There is no grid point right at the 
interface; the difference in material my be accounted for by 


introducing an interface diffusivity 


4 = (4.20) 


in which k,denotes 2 harmonic mean conductivity (46) 


material 2 


Figure 4.1 


interface material 1 
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Interface properties 
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(44) 


and 956, 


° 
Sip. * Pet 
9c, , =e eae (4.12) 


With the use of these quantities, the continuity of heat fluxes at 
the interface may be satisfied implicitly. 
4.1.3 kinematically Consistent Velocity Fields 

It has been well established that finite difference 
approximations for conservation of mass, vorticity, and energy are 
relatively easy to derive for a system in Cartesian coordinates. In 
cylindrical coordinates, however, they are more difficult to derive 
because the transport velocities (normal velocity components on grid 
volunes), if improperly expressed in terms of stream functions, my 
lead to a situation in which the conservation of volume flow is not 
satisfied on the grid volume. 

Parmentier and Torrance [48] derived kinematically consistent 
velocities for use in cylindrical coordinates. The derivations that 
follow are credited to their work but modified slightly in order to 
account for the nonuniform grid system used in the present study. 

It is imperative that kinematically consistent velocities must 


satisfy the circulation theorem both locally and globally 


{sds =f WaT (4.13) 
1 
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In this equation s and 1 represent respectively the area and the 
perimeter of the region of integration. In a dimensionless form, 


equation (4.13) can be expressed as 


fader =A f udz +f var (4.44) 
s 

where A is aspect ratio, H/R;. Approximation of the integrals in the 
above equation leads to the following expression for the area- 


averaged vorticity in the grid volume: 


AUP jen 7 Uf Mey eas vb gl 
potest Pr a tad M4-1,5°%9 
og oe (4.15) 


ce) 


Consider an interior grid point in a nonuniform grid system 
shown in Figure 4.2. Since the faces of the grid volume are located 
midway between the mesh points in each direction, approximations for 
tangential velocities in terms of stream function my be obtained 


from equations (3.20) and (3.21) as 


(4.16) 


(a3) 


Figure 4.2 Notations for an interior grid volume in a nonuniform 
cylindrical grid system 
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Velocities at grid points may then be expressed in terms of these 


velocities as 


Caer 
Oe oe is iin 
ge 
(ary 
© git 
a os tag aisig 
se) 
Upon substitution of equation (4.16), equation (4.15) becones 
i 
(4.18) 


This 15 the finite-difference form of the stream function equation 
(3.19). Here the circulation theorem has been used for deri vation; 
Taylor's series my be used to derive the same equation, as is not 
unexpected. 

For an incompressible medium, the normal velocity components 


(transport velocities) must be defined with the following relations 


(4.19) 
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50 that the volume flow is conserved on the grid volume. The 
asterisked stream functions refer to the stream functions at the 
corners of the grid volume; see Figure 4.2. These corner stream 
functions are related to the stream functions at grid points by the 


relation 


* Yi,5* Yin ed * Mis je 
Wg 7 age) + 9 (Aa sey (4.20) 


In order for the mass flow to satisfy the conservation principle, the 
weighting functions, aj and 2j, must take on those values listed in 
Table 4.2. 

Combining equations (4.19) and (4.20) gives expressions for the 


‘transport velocities as 


A 
Rig 

(4.21) 
a =e Pitty j Yin, Viel, jot ~ Yi-2, fol. 
Ys Teta Jy + g(t Pe eh 


On the other hand, the grid-point velocities may be derived by 


combining equations (4.16) and (4,17) as 


Table 4.2 Weighting coefficients in equation (4.20) 


J ay i 
1 a/4 Va 
2 TAZ S/2 
3 1/20 9/20 
4 15/28 | 13/28 
5 19/36 | 17/36 
6 23/44 21/44 
7 27/52 25/52 

28 ve ve 
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(4.22) 


The derivations given above permit the vorticity and velocity 


boundary conditions to be written as follows: 


Vorticity Boundary Condi tions: 


2-0 


2H 
“i9 re 


ar jaz 


(4.23) 


(4.24) 


(4.25) 


(4.26) 


Note that the three vorticity equations given above are first-order 
accurate. However, they are the only ones that are consistent with 
the interior vorticity field and are able to satisfy the circulation 
theorem given earlier. 

In the numerical computation that will be made later, equation 
(4.22) will be used to compute the grid-point velocities. Equation 
(4.21) will be used to compute the transport velocities for use in 
‘the convection terms in the vorticity transport equation (3.18) and 
the energy equation (3.22). 

4.1.4 Generation of Nonuniform Grid System 

For the convection problem under investigation, it is expected 
that there will be a boundary layer near the wall anda recirculating 
Flow in the core region. It is desirable to use a nonuniform grid 
system such that a denser grid is formed in regions where both the 
Flow and temperature fields change more rapidly. The method used to 
create a nonuniform grid in this work is the coordinate stretching 
technique. 

Two stretching functions are used in this work. In the z 


direction, the stretching function used is (49) 


n(z) = 0,5{1 + tanl(2z-1) F8I/tan $5} (4.27) 


This function enables a mapping of z to ny © being a deformation 
parameter--its presence controls the spacing of grid lines. 
In practice, the velocity and temperature fields are solved 


nunericallys it will be convenient if the indices used to represent 


a 


positions in the numerical computation can be related to the actual 
positions in the physical domain. The inverse of equation (4.27) is 
thus more useful and given as 


2(n) = 0.5 +45 tan 


C(2n-1) tan $4) (4.28) 


Examples of coordinate transformation constructed for ¢ = 0.2 and 0.9 
are shown in Figure 4.3(a). The transformation of e = 0.2 is close 
to being linear in that a uniform grid size in n is transformed to a 
uniform grid in 2, For the present study, use of © = 0.2 is 
adequate. 


A different function is used to stretch grids in the r direction 


gir) = ale” - 1) (4.29) 


This time a and b are deformation parameters, and for the present 
work, 2 and b are taken to be 0,044 and 3.1666, respectively. These 
values generate a grid system that has a coarse grid near r= 0. The 
grid lines become denser as they move close to the wall, see Figure 
4.3(d). 


The inverse of equation (4.30) s 
PAs ee 
rehinete (4,30) 


This equation enables the positions of r be found if ¢ values are 


given, 
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ching functions in (a) z direction and 


s 
(3) r direction 


Figure 4.3 
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A simu] taneous use of equations (4.28) and (4.30) creates a 
41X21 grid net as shown in Figure 4.4, This is the grid system used 


for numerical computation in the present study. 


4.2 Solution Strategy 

It is necessary to make decisions on (1) how to discretize the 
partial-differential governing equations given in Chapter III, and 
(2) how to solve these discretized equations numerically. 

The first decision, in fact, has already been made in Section 
4.1, As has been discussed there, a control-volume approach is able 
to yield a kinematically consistent velocity field. In addition, in 
formulating the boundary conditions for meeting continuous 
‘temperature and heat-flux conditions at the wall/insulation 
interface, an interface-diffusivity concept, equation (4.10), mas 
also been employed. This concept is again developed based on the use 
of control-volume approach. This control-volume approach will be 
used consistently in formulating other finite-difference equations in 
this work [46,50]. 

As for the second decision on the selection of solution methods, 
‘two methods will be used for solving two kinds of equations. For 
example, the vorticity transport equation (3.18), the water energy 
equation (3.22), and the wall energy equation (3.25) a1? belong to 
‘the same types namely, they are parabolic in time but elliptic in 
space. Such equations can be solved with several methods; the one 
that is chosen here is the Alternating Direction Implicit (ADI) 
method of Peaceman and Rachfora [51]. In this method, each time step 


is split into two half steps; the space derivatives are approximated 


Figure 4.4 Grid system in water used in the present study 
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with both implicit- and explicit-difference schemes, alternating in 
each direction. Then the finite-difference equations can be 
expressed ina system of linear equations which can be organized in 
‘tridiagonal matrix forms suitable for solution with an efficient 
Thomas algorithm [50]. The ADI method is superior because it allows 
the use of large time steps to solve problems with a large number of 
equations (521. 

It must be noted, however, that the water energy equation and 
the vorticity transport equation are both algebraically nonlinear in 
‘the unknowns, and this is because of the presence of convection terms 
which contain these unknowns in the coefficients expressed at new 
time. Linearization of these equations is necessary, and the method 
chosen is to evaluate these coefficients at old time--commonly known 
a5 lagging the coefficients [50]. In practice, the convection terms 
in the energy and vorticity transport equations are evaluated at old 
time, whereas the diffusion terms are split following the ADI method 
mentioned earlier. This hybrid scheme makes it possible to resolve 
the difficulty of nonlinearity while maintaining the transportive 
nature in those terms characterizing heat diffusion [36]. 

The stream function equation (3.19) is of elliptic type. For 
‘this equation, methods for solution are different and usually fall 
into two categories: direct and iterative. The direct method is 
limited to the solution of systems of small numbers of equations; the 
large round-off error inherent in this method also makes it 


unfavorable. The iterative method is therefore chosen. In addition, 


a successive-over-relaxation (SOR) is employed in order to speed up 


convergence in the iteration. 


CHAPTER V 
FINITE DIFFERENCE EQUATIONS FOR GOVERNING EQUATLONS 
‘AND BOUNDARY CONDITIONS 


The decisions made in Chapter IV on solution strategies will be 
used as guidance here in developing finite-difference approximations 


of the governing equations and boundary conditions given in Chapter 


MI, 
5.1 Fini te-Difference Formulation for 
nergy Equation in Wall 
The energy equation for wall has been given in Chapter III as 
follow 


2, 
alat,) — a°(at,) 
Lae B 
+ (3.26) 
ae 
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The associated boundary conditions in the wall are 


(3.27) 
ret, Ty(O,2,t) = Tylzt y 


a7 
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By integrating equation (3.26) over the control volume about 


point (i,j) (Figure 5.1) for a half-time step 


At 
nw te aT, nw tho a(aT.) 
z 3 o 2 
ss Cre), gatdr gaz = IS laa, jdtdr gaz 


at 
nw te St acar.y 

os ie 
set 


nw te St. a(t) 
AIT Baa, jdtdrgdz (5.1) 


it is possible to derive the finite-difference approximation of the 


energy equation as 


(5.2) 
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Figure 5.1 Notation for control-volume method 
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For this half-time step, the temperature is expressed explicity in 
the r direction but implicitly in the z direction. According to the 
ADI method descrived in Chapter IV, for the next half-time step, the 
expression for these temperatures will be switched in that the 
temperature Will be expressed explicitly in the z direction but 
implicitly in the r direction, Further development of fini te- 
difference formulation in the wall will depend on how the wall is 
constructed. 
5.1.1 Bare Wall 

For a bare wall without insulation on it, the dimensionless 
diffusivity, a, that is the ratio of wall-diffusivity to water- 
diffusivity, is constant throughout the wall. Then, the tempera ture 
at the interface of two neighboring control volumes (Figure 5.1) may 
be approximated by an arithmetic mean 

1 


gd BT ga Tg) (5.3) 


and the temperature gradients approximated as 


(5.4) 


With the substitution of these relations and dropping subscript s of 


T for simplicity sake, equation (5.2) is changed to 
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1 1 1 
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This equation can be expressed compactly as 


1 1 1 
wed wed keh 
2 eye 2, z 
(my Tags + Ret yy ty IT SSP + Cae IT Ly 
k 2 k k 
= MCyg ~ CygITE jaa + em AygITH 5 + (yg + CyyITE yy (568) 


where coefficients C); through Cjq are defined in Table 5.1. 

In the derivation given above, equation (5,6) was obtained by 
integrating equation (3.26); this same equation (5.6) may be derived 
by using Taylor's series, as is not unexpected. 

For the second half-time step, the finite-difference equation 
for T can be derived as 


kt k+l k 
Chg * Cag jan * OBE Ag ITE + (05 - Sw 5a 
1 il et 

keh ic wet 


= (004 IT; (5.7) 
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Table 5.1 Definitions for C coefficients 
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Table 5.i-continued. 
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where ay =u, = [up] 


Yq and v, are defined in equation (4.21) 


rj, ar; are defined below and shown in Figure 4.2 


other variables are defined in NOMENCLATURE 


Wy (24g > 24/2 
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5.1.2 Insulated Yall 

For a tank wall with insulation on it, temperatures for interior 
grid points in wall (or insulation) will be treated just like those 
grid points in a bare wall but will use correct properties to account 
for material variation. For those points adjacent to the 
wall/insulation interface, however, special care must be taken to 
account for material property variation across the interface. 

Consider the grid point located to the right of the interface 


shown in Figure 5.2(a). It is obvious from equation (5.2) that 


(5.8) 


Here aj and aj1) stand for dimensionless diffusivity of material 1 
because the dashed-line boundaries separating the control volumes to 
the right of the interface are a11 located in material 1. The 
diffusivity aj41p. however, is different; this diffusivity is 
defined as the ratio of the ‘interface’ diffusivity (equation (4.10)) 
to water diffusivity. 

The same distinction must be made in temperature expressions. 


White Ty, 5-1 cam be approximated as the arithnetic mean of 
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Figure 5.2 Treatnent of grid points adjacent to interface 
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neighboring temperatures Ty; and Tj, 5.1 Tj, je must be expressed 
differently. In order to satisfy continous heat fluxes at wall/water 


interface, the ‘interface’ temperature T, Must meet similar 


inde le 
requirement of flux continuity, thus satisfying 


ee (5.9) 
Ley rr +e 
1 


With the substitution of these new relations, equation (5.2) takes 


the following form: 


wed wed ed 
2 
Cay Tangy GEE Cay # mya ITG, FP + Cobra T5615 


see ae 
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peewee 
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where Cy, through Cy4 are again defined as in Table 5.1. Those CT, 
through C¥s are defined as foltows: 


* 2 o°B 
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13" ie 
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In these relations, a” is the interface diffusivity ratio, as has 
deen defined previously. 
For the second half-tine step, the finite-difference equation 


for T can be derived as 


hal 2 ke ke 
(CastaalT ja * Ge Maat IatCre tials HCI Chg TE a 
ae et ed 
2 aby Tyas + Ge tay tel yay P+ lata yy (5.12) 


The derivation given above for an insulated wall is for a grid 
Point located to the right of the interface (Figure 5.2(a}). For the 
grid point located to the left of the interface, Figure 5.2(b), a 

slight modification of the formulation is in order. This time aj. Ip 
strides across two different materials and thus becomes the interface 
diffusivity. Correspondingly, Tj, 5.1), 18 the interface temperature 


25 defined previously. For this grid point, the first half-time step 


and second half-time step fini te-difference equations are 


respectively 
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wh wd 
2 2 
+ Ge ee ial, 3 - (ayaa gy 


Ps 
= (00 IT; (5.18) 


where the definitions for cosfficients remain unchanged. 

The finite-difference formulation for governing equations in the 
wall is now complete. Attention is now directed to the boundary 
conditions in the wall. 

For grid points located at the air/wall interface, Figure 5.3, 
‘the boundary condition for T is 


ig ho Mets) 


For grid points along the top and bottom boundaries of the wall, the 
adiabatic condition there requires that the following condition be 


met 


(5.16) 


Along the wall/water interface, the boundary temperature must satisfy 


Ti - Tinterface Set 


These boundary conditions may be incorporated into the finite- 


difference energy equations derived earlier for the wall and to 


Tinterface 
air/watl wall/water 
interface ‘interface 


Figure 5.3 Wall region temperature boundary conditions 
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Simplify those equations that contain grid-point temperatures 
adjacent to the boundaries. These difference energy equations are 
‘then solved and once the temperatures at interior points are found 
numerically, the boundary temperatures may be recomputed using the 
boundary conditions given above. 
5.1.3 Sumary in Matrix Forms 

The finite-difference governing equations derived above may be 
expressed in matrix forms, as will now be shown. In general, the 


‘temperatures in the wall are found by solving the matrix equation 


ATG (5.18) 


Matrix A and vectors 8 and C in this equation vary with time steps 
and also the way the wall is constructed. 


(1) For a bare wall without insulation on it, at the first half 


time step, 
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(5.190) 


(5.19¢) 


At the second half-time step, 
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where C)) through C14 have been defined in Table 5.1. It is clear 
that both A matrices given in equations (5.193) and (5.20a) are 
‘tridiagonal, a form enabling the matrix equation (5.18) to be solved 
by efficient Thomas algorithm, as mentioned in Chapter IV. 

(2) For an insulated wall, at first half-time step, matrix A is 
unchanged from that of the bare wall at the same time-step. 
Likewise, vector 8 follows much the same formas before, except for 
‘those two grid lines adjacent to the wall/insulation interface. Say 


the interface is located between j = j-l and j = j+l, then 


65 


* oak stotaet yk 
31s aga, 5 a3 a5 MT, jar 


k 2 
(yg Cag Ta, jan HE 


=6 sik ES a nok yk + get yh 
asian, jae Cas Cis Cag aT, MSI MT, jet 
(5.2la) 
isan 
# ot qk 2 + Ore ae k 
(rsa lTs, joa MEE Cas Cag Cas ITG, 5 Mas 1g ITY, jan 


gt yk 2 
(311, 5a ae 


Jig cet ih k 
13761314 Cis Tian, saga, en 
(5,21b) 


At second half-time step, vector 8 fs unchanged from that of the bare 


wall at the same time-step. Matrix A is changed as follows: 
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(5.22) 


5.2 Finite-vifference Formulation for Energy Equation in Water 


The energy equation in water has been given in Chapter III as 
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The associated boundary conditions in water are 


(3.23) 


fs 
aero 


Tp(L,24t) = Tylzst) 


The Finite-difference energy equation in water my be derived 


following the same procedure used before for the wall and the result 
is 


(5.23) 


This equation is derived for temperatures of water at the first half- 
time step; here the temperature is expressed explicitly in the r 


direction but implicitly in the z direction. Velocity coefficients, 


being taken at old times, are ‘lagged’ in the computation, as 
described in Chapter IV. 

A Second-upwind differencing scheme is used to treat the 
convection terms in equation (5.23). Referring to Figure 5.4, 


1 

a) > uT) = PUluy UIT 4g 5 HUgtlug IT). 9 
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Here Us, Uns Ve and vy refer to old-time transport velocities defined 
by equation (4.21). With the use of these two relations, equation 


(5.23) is changed to 
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where C11, Cy, and C\¢ through C,, have been defined in Table 5.1. 
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Figure 5.4 
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Second upwind difference scheme 
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For the second half-tine step, 
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The formulation of the finite-difference energy equation for 
water is now complete. Attention is now directed to the boundary 
conditions in water as follows: 


For grid points at the centerline of the tank, Figure 5.5, 


a (5.27) 
For grid points located at the top and bottom boundaries, 
Teh Tas 
(5.28) 
Tha Teas 
For grid points at the wall/water interface, 
Tia * Tinterface (2) 


Once again, these conditions are incorporated into the fini te- 
difference energy equations for water, and these equations are solved 


for interior temperatures. Once these interior temperatures are 
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Tinterface 


wall/water 
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Figure 5.5 Water region temperature boundary conditions 


2 


found, the boundary temperatures of water can be obtained by using 
these boundary conditions given above. 

{tis again desirable to sumarize the governing equations 
derived above for water in matrix expressions. In general, the 


interior temperatures of water are found by solving 


Re 


8 (5.30) 


The matrix A and vectors 8 and C are defined differently for 


different time steps. For example, at the first half-time step, 
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Coefficients C11, Cy2 and Cyg through Co) have been defined in Table 


5.1. Again note that the bio A matrices given are tridiagonal. 


5.3. Iterative Solution for Two Energy Equations 


Soth energy equations for wall and water must be solved 
simultaneously to satisfy continuity of temperatures and heat fluxes 
at the wall/water interface. It is anticipated that many iterations 
are needed to obtain convergent temperatures at each time step. In 
this effort, the already convergent, previous tine~step tempera tures 
(for example, temperature vectors in equations (5,19) and (5.31b)), 


are assigned fixed at each iteration and new temperatures at k+(1/2) 
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and k+l at the same grid locations are computed, iterated, and tested 
until they converge. Each time during iteration, the interface 
temperature is modified by using equation (4.9) and the temperature 


convergence is hinged on a convergent heat flux 


(5.33) 


where the interface grid locations are indexed as i, i=l through I. 
Heat fluxes gp and qg refer to those in the water and the wall, 
respectively. 
5.4 Finite-Difference Formulation for 
jor ticity Transport Equa tfon 
The vorticity transport equation has been given in Chapter III 


as 


aeplve) + Spf) = 2B TE a PrtA seth Pdr gu)) 
apes 
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Its finite-difference form for the first half-time step is 
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where the vorticity is expressed explicitly in the r direction but 
‘implicitly in the z direction. For the second half-time step, the 
vorticities are switched between their explicit and implicit 
expressions, giving 


Mey a ke kel 
(ean eh gen * Cae Seales + (Cag) ha 


wed wed 
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Fe Hada gt bagdogan 2s ~ Cg( THE, = THE) 
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where Cog through Cyq have been defined in Table 5.1. 
Since the driving force for buoyancy-driven flow in water is 
induced by the water temperature variation, the water energy equation 

‘must be solved prior to the solution of the vorticity transport. 
equation. The temperatures of water solved from the energy equation 
are used in the vorticity transport equation, and these temperatures 
are kept constant for eacn ‘full’ time step in the iteration. 

The vorticity boundary conditions have been given in Chapter IV 


as 


(4.23) 
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Sra (4.24) 


TPM ded 


yt ats (4.25) 
oe 

The boundary grid-point locations are shown in Figure 5.6. The three 

boundary conditions given above are for solid boundaries; at the 

centerline, the vorticity is zero because of equation (3.14), in 

which the terms Vj. and 3V,/ar both vanish. 

In practice, with the exception of the known vorticity along the 
centerline, the vorticity values along the solid boundaries are 
unknown a priori; an iteration must therefore be used to find them. 
The approach used in this dissertation is (i) make a guess of tne 
boundary vorticities, (ii) solve vorticities at all interior points 
in water using the vorticity transport equations (5.34) and (5.35), 
(iti) use the vorticities found in step (ii) in the stream function 
equation (3.19) to obtain new stream function values, and (iv) use 
these new stream function values in equations (4.23) to (4.25) to 
compute a refined set of boundary vorticities. 

In order to safeguard the computed vortici ties from converging 


to erroneous 


Its, as warned by Israeli [53], the boundary 


vorticities computed using equations (4.23) to (4.25) are fine-tuned 


top 
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Figure 5.6 Water region vorticity boundary conditions 
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80 
using the following formula: 


=o) +) (5.36) 


Im wits equation, o) refers to the boundary vortici ties conputed from 
equations (4.23) to (4.25); ¥, stands for a finite-difference approxima- 
tion to the normal derivative of the stream function evaluated at the 
boundary. The quantity « represents an iteration parameter whose value 
is taken to be 0.5 in the present study. The superscripts 1 and 1#1 are 
sed to designate iteration index. Equation (5.36) enables computation 
of convergent boundary vorticities while inducing the no-slip condi tion 
‘that is necessary for any solid boundary. 

In practice, the boundary vorticities at each half-time step are 
obtained by interpolating the old and new vorticity values as follows 
[54]: 


= (a) + ally (5.37) 
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were ol! refers w the vorticity caleuated fron equttion (5.38). 
The convergence of the vortiities conputed at each grid point is 


checked by using [55] 


Max 
boundary 


(5.38) 


This convergence criterion will be reexamined in Chapter VII. 
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When the vorticity boundary conditions are incorporated into the 
Finite-difference equations (5.34) and (5.35), the resulting 
equations can be expressed ina matrix form 

Re = 65 (5.39) 
where matrix A and vectors B and C are defined differently for 


different time steps. For the first mlf-time step, 
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At the second half-time step, 
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where Coq through Cyg have been defined in Table 5.1. Once again, 


both A matrices are tridiagonal. 


5.5 Finite-Difference Formulation for Stream Function Equation 


The stream function equation has been given in Chapter III as 
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ir; 


Ms 


y+ Natio se (3.19) 
Since the solid boundaries and the centerline of the tank can be 
considered as one of the streamlines, a stream function value zero is 
chosen for this streamline, As has been mentioned in Chapter IV, 
‘this stream function equation (3.19) will be solved by using a 
successive-over-relaxa tion (SOR) method, as wil1 now be shown. 


The finite-difference approximation of equation (3.19) is 
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where Cg through yp have been defined in Table 5.1. The Gauss- 


Seidel iteration is used to write 
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where the superscript 1 refers to the iteration index. The computed 
stream function value is over relaxed by using the formula 


si 


WTS = = eho aatet (5.44) 


ar) 
where X represents the relaxation factor and a value of 1.7 is tested 
for use in the present study. In the numerical computation carried 
out in this work, the stream function values at each grid point are 
computed repeatedly until these values converge and satisfy the 


following criterion 


vax | ted | ¢ 04 (5.45) 


5.6 Solution Algorithm 
To close this chapter, a solution algori thn is provided as 


follows: 


a) 


(2) 


a) 


a) 


(s) 


(6) 


a) 


(a) 


Use the temperature of the previous time step as an initial 
guess. (For the first tine step, the initial tenperature is 
Used to initiate computation.) 

Compute a new set of temperatures at each grid point using 
equations (5.18) and (5.30). 

Check heat-flux continuity at wall/water interface by using 
equation (5.33). If this condition is not met, then modify 
the interface temperatures using equation (4.9) and return 
to step (2) above. 

Use the vorticities of the previous tine step as an initial 
guess. (For the first tine step, the vorticities are 
assumed to be zero everywhere to initiate computation.) 
Solve the vorticity transport equation (5.39) at interfor 
grid points. 

Solve the stream function equation (5.46) using the new 
vorticities obtained in step (5) above. 

Compute new boundary vorticities using stream functions 
found in step (6) above. Check the convergence of boundary 
vorticities using equation (5.38). If the convergence 
criterion is not met, then modify the boundary vorticities 
using equation (5.36) and return to step (5) above. 
Calculate grid-point velocities using equations (4.22) and 
(4.26). 


The procedure given above is for one (full) time step. The 


whole procedure is repeated for the next time step until the desired 


time is reached, 


CHAPTER VI 
RESULTS AND DISCUSSION 


The enclosure under investigation has an aspect ratio 4; it 
measures 60 cm high and 15 cm radius. The enclosure wall is made up 
of 12 carbon steel and insulated with soft rubber, both being 0.4 cm 
‘thick. Thermophysical properties of these materials are given 
together with water in Appendix A. The initial (dimensionless) 
temperature distribution in water is given as shown in Figure 3.1; 
the initial temperature of the wall is unity throughout. The actual 
temperatures corresponding to the dimensionless temperatures 0.2 and 
1,0 marked in the figure are 20°C and 60°C, respectively. The 
ambient air temperature is kept constant at 10°C. 

In the numerical computation, a grid system of 41X21 points is 
sed for the water. Grid systems of 41X10 or 41X6 are used 
respectively for the wall with, or without, insulation. A (full) 
time step used in the computation is 2 seconds. Selection of the 
grid size and time step is subjected to tests later (Chapter VIL). 
The numerical results presented in this chapter were obtained using a 
Harris 800 computer and computed with double precision. 

In order to view the flow and temperature fields in water, 
isotherms are plotted together with velocities and stream functions; 
see, for example, Figure 6.1. The spaces provided for velocities 


(middle figure in each set) are too small to plot vectors; only line 
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segments are thus shown, With the help of the stream function plots 
given on the right, the direction of the velocity can be 

identified. Specifically, a positive stream function value 
represents a counterclockwise circulation loop and vice versa. Soth 
the velocities and stream functions plotted are normalized by their 
instantaneous maximum values found in the flow field. The thickness 
of the wall in the isotherm plot (left figure) is tripled in order to 


give a clear view of the isotherms. 


6.1 Flow and Temperature Fields in Water for 
Dpifterane bestgne oF te WaT 

Prior to the presentation of nunerical results for comparing the 
Flow and temperature fields resulting from wall design variations, 
physical mechanisms governing the flow and heat transfer are examined 
and given below. 

(4) At the Tower section of the tank where the wall temperature 
is initially higher than water, heat always transfers from the wall 
to the water during the early stage of cooling. At this time, since 
the convective resistance at the water side is always smaller than 
‘that at the air side, more heat tends to flow to the water from the 
wall. 

(ii) Since the ambient air fs always at a lower temperature, 
the wall will lose heat to the air once cooling starts. dn the other 
and, for the tanks under investigation, the water at the top section 
of the tank is always at a temperature higher than that at the 
bottom; heat will thus leak from the top to the bottom via the tank 


wall. This heat-leak problem is particularly pronounced for a tank 
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wall made up of @ good thermal conductor. For sucha wall, the 
Physical mechanism Teading to the heat leak is that the wall, being 
heated by the not water above, is able to maintain a temperature that 
is slightly higher than the water at the same level; a convective 
Joop thus forms at this level. As to how far this heat will 
Penetrate downward along the wall depends greatly on the 
effectiveness of heat transfer in the wall as well as the amount of 
heat that is simultaneously lost by convection from the wall to the 
outside air, The difference between these two controls the heat 
penetration length. 

6.1.1 Uninsulated Tank 

For the ease of reference later, this case of tank wall without 
insulation on it will be called Ni for short. For this case, because 
of the physical mechanism (1) described above, clockwise circulation 
loops form in the lower section of the tank; see Figure 6.1(a). The 
wall being highly conductive, the temperature gradient across the 
wall is negligibly smi1, a situation contributing to the horizontal 
isotherms viewed in this region. 

For this NI case at 180 seconds, the temperature of the lower 
wall drops toa level where the lower circulation loops in water 
Start to weaken. In the mean time, the top wall cools down due to 
heat loss to the ambient air; circulation loops thus form in this top 
region. These top and bottom loops move in opposite direction as 


revealed by the stream function values marked. 


(2) 60 seconds (8) 100 seconds 


Figure 6.1 Isothora, velocity Field, and strenaline for tank without insulation 
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{e) 7200 seconds (1) 10000 seconds 


igure 6.1-continued. 
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A short while later at 300 seconds, tne lower circulation loops 
nearly disappear; the top circulation loops, however, still 
Persist. In fact, because of the continuous heat loss from the wall, 
‘the upper circulation loops intensify and are now able to penetrate 
deeper and deeper into the lower region; see Figures 6.1(d), (e), and 
(f). A quasi-steady state is finally reached, and the stratification 
of water in the tank is also destroyed as shown by the isotherm 
locations plotted; also see the centerline water temperatures plotted 
in Figure 6.2. 
6.1.2 Tank Insulated from the Outside 

Placing insulation over the outside of the wall (referred to 
later as EI case for short) increases the initial heat transfer from 
the wall to water. As expected, the heat loss from the wall to the 
ambient air is now reduced; as a result, the upper circulation loops 
cannot develop as rapidly as in the previous NI case compare Figures 
6.1(b) and 6.3(b). The heat leak via the conductive wall, however, 
becomes more serious; see the circulation loops plotted in the lower 
section of the tank at 300 seconds, Figure 6.3(c). A comparison of 
the isotherm locations also indicates that the stratification in 
water decays more rapidly in this case. The centerline water 


‘temperatures given in Figure 6.4 also serve to illustrate this point. 
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Figure 6.2 Centerline water temperature variation with time for tank 
without insulation 
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Figure 6.3-continued. 
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Figure 6.4 Centerline water temperature variation with time for tank 
with exterior insulation 


99 


6.1.3 Tank Insulated at the Inside 

For a tank wall insulated at the inside (II case), the heat leak 
from the hot water to the cold water is inhibited because of the 
insulation. This reduction of heat leak results in smaller 
circulation loops both at the top and bottom sections of the tanks 
see Figure 6.5. At very long time (3,600 seconds), small 
counterclockwise loops also develop at the bottom section of the 
tank, a result of the heat loss from the lower wall to the ambient 
air, Stratification in water is well preserved in this case, as 
demonstrated by the isotherms plotted in Figure 6.5 as well as the 


centerline water temperatures shown in Figure 6.6 


6.2 Heat Fluxes at Nall/Water Interface 
An attempt is also made to calculate the heat flux along the 
wall/water interface. Fourier's law, in finite-difference form, is 
used to evaluate the heat flux at the water side, and the results are 
Plotted in Figures 6.7 through 6.9, In these figures, a positive 
heat flux value stands for a heat flow from water to wall and vice 
versa, The z coordinate is dimensionless, while the heat flux is 


dimensioned, 
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Figure 6.6 Centerline water tenperature variation with time for tank 
with interior insulation 
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6.2.1 Heat Flux Distribution for Ni Case 

The observations made in section 6.1 relating to velocity and 
temperature fields can be supported by the heat fluxes presented in 
this section. Take the NI case at 60 seconds (Fig. 6.7) for 
example. Heat flows from the wall to water in the lower section of 
‘the tanks in the top section, the heat flux direction is reversed. 
For longer times at 180 and 300 seconds, the lower heat fluxes 
diminish rapidly. There are points at about z = 0.2 and 0.6 where 
‘the heat fluxes reach their maximum values. This may be alluded to 
the heat penetration length discussed in section 6.1. The heat 
‘transfer in the tank wall is complicated because, on one hand, the 
wall receives heat from the hot water above, while on the other hand, 
the wall also dissipates the heat that it receives to the ambient 
air. The problem is further complicated by the fact that, at the 
water side, the temperature is stratified. 

For times past 300 seconds, heat loss from the water comes 
mainly from the top. The heat flux curves plotted in Figure 6.7(b) 
Serve well to illustrate this point. 

6.2.2 Heat Flux Distribution for El Case 

For a tank wall insulated from the outside, the heat transfer 
from the lower wall to water is increased; note the scale change for 
9 in Figure 6,8. The heat loss from the top section of the wall is 


reduced because of the insulation. 
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Figure 6.7 Heat flux distribution for tank without insulation 
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Figure 6.8 Heat flux distribution for tank with exterior insulation 
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6.2.3 Heat Flux Distribution for {1 Case 

For a tank wall insulated at the inside, the heat loss from the 
top section of the tank is the lowest of the three cases studied; see 
Figure 6.9. 

To facilitate comparison of the heat fluxes computed for the 
three cases above, a composite diagram is constructed as shown in 
Figure 6.10. A change in scale for q should be noted for Figures 
6,10(b) and (c). 


6.3 Convective Coefficients at the Air Side 

The convective heat loss from the wall at the air side can be 
evaluated by using equation (4.5). Newton's cooling law can be used 
together with this equation to find the local convective coefficient 
28 plotted in Figure 6.11. It is difficult to precisely pinpoint why 
the convective coefficient varies as plotted because, from equation 
(4.5), the local heat flux from a nonisothermal wall is a complicated 
function of the local temperature excess over the ambient 
‘temperature, and the gradient of this temperature excess along the 
wall, among others. Some basic features can nevertheless be found 
from these curves as described oelow. 

For all the three cases studied, the convective coefficient is 
always the largest near the bottom end of the wall, a result of the 
smal] boundary-layer thickness there. This coefficient drops off 
rapidly as it moves up along the wall. The variation of the 
convective coefficient is somewhat irregular at higher elevations; 


however, the magnitude of this coefficient is still reasonable, 
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Figure 6.9 Heat flux distribution for tank with interior insulation 
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gure 6.10 Heat flux distribution for three designs of the wall at 
different times 
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Figure 6.10-continued. 


uz 


1.9) 
NO INSULATION 
2 ae 
300s 
#3600 5 
5p 
ool + * - 
10 
ob 
Py cena figs 2 
INTERIOR INSUL. 
z 5 
300s 
3600 5 
o.5b 
1 coe 


untae 


Figure 6.11 Heat transfer coefficients for three designs of the wall 
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For all the cases tested, the convective coefficient varies from 3.5 


to 5 wint*c, 


6.4 Velocity and Temperature Distributions 
ete tomer ego 


The velocity and temperature profiles in the boundary region at 
‘time equal to 60 seconds are plotted in Figures 6.12 and 6.13. The 
elevation z and the radial position r in these plots are 
dimensionless, whereas the axial velocity is dimensioned. The 
dimensionless temperature plotted is measured in excess of that 
temperature measured at the centerline and at the same elevation. 
Some trends can be seen from these curves. 

As shown in Figure 6,12, of the three cases studied, the EI case 
always has the largest velocity in the boundary layer; the velocity 
for the II case is the smallest. Similarly, in the temperature plots 
given in Figure 6.13, the dimensionless water temperature is always 
‘the largest for the EI case and the smallest for the II case. 
Because of the short time (60 seconds) chosen for these plots, the 
lower half of the wall is always at a temperature higher than water, 
‘the velocities of water in the boundary region are thus positive 
close to the wall. 

It is informative to examine the magnitude of the water velocity 
‘in the boundary region. The maximum axial velocity of water 
calculated for EI case is 0.38 cm/s. Also, the thickness of the 
velocity boundary layer where the flow velocity is positive accounts 
for about 5 to 10% of the tank radius. These magnitudes appear to be 


in reasonable agreement with the 


perimental data reported by 
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Figure 6.12 Velocity profiles at various elevations at 60 seconds 
for three designs of the wall 
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Figure 6.13 Temperature profiles at various elevations at 60 seconds 
for three designs of the wall 
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Figure 6.13-continued, 
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fies and Miller [18] who used a laser Doppler velocimeter to measure 
the water velocity near the wall in the cold region of a 
thermocline. This region corresponds to the lower section of the 
tank in the present study; use of their data for comparison in Table 
6.1 {5 thus justified, 

A close examination of the velocity plots given in Figure 6.12 
reveals that, for tank elevations above about 8/40, the velocity 
profiles consist of tio branches: an upward-pointed velocity region 
close to the wall, and a downward-pointed velocity region away fron 
the wall. The tanperature profiles given in Figure 6.13 can be 
divided in a similar fashion; the dimensionless tenperature close to 
the wall is positive, while this temperature away from the wall is 
negative. These peculiar velocity and temperature distributions are 
commonly referred to as flow reversal and temperature defect, 


respectively, as will now be discussed. 


6.5 Flow Reversal and Temperature Defect 
A detailed study of the flow reversal and temperature defect 


requires a close examination of the flow field in the lover section 


of the tank. Take the NI case for example; an enlarged view of the 
flow field near the wall up to an elevation of grid line 21 is given 
in Figure 6.14. The directions of the line segments, representing 
velocities, are slowly turning in the region away from the wall. To 
support this observation, the velocities at two elevations (i = 13 
and 11) are replotted as shown in Figure 6.15, in which each velocity 


vector is decomposed into its z and r components. Ata lower 


ug 


Table 6.1 Comparison between numerical and experimental data 


Expertimenc* Present study** 


Axial velocity inside | 0.45 cm/s (maximum) | 0.38 ca/s (maximun)} 
‘the boundary layer 


0.25 cm/s (typical) | 0.02 ~ 0.38 ca/s 


Core axial velocity | 0.01 cm/s (typical) | 0.005 ~ 0.056 cm/s 


Boundary layer 5 ~ 10% of radius 5 ~ 10% of radius 
‘thickness 


i 
\ 
f 


data (18] given are for observations made within the first 3 
minutes of the experiment for EI case 


Ie EI case at 60 seconds 
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Figure 6.14 An enlarged view of the flow field in the lower section 
of the enclosure in Fig. 6.1(a) 
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Figure 6.15 Decomposition of velocity vectors into their z and r 
components 


ize 


elevation of i = 11, the r-component velocities increase with 
decreasing r (moving away from the wall), contributing to a rapid 
rotation of the velocities at this elevation. This velocity rotation 
and eventual reversal are largely responsible for the wedge-1ike 
Streamlines viewed at point A in Figure 6.1(a). As expected, the 
radial component of the velocity causes flow recirculation; the 2- 
component velocity gives rise toa flow reversal. 

Study of flow reversal and temperature defect is facilitated by 
simu] taneously viewing the velocity and temperature profiles plotted 
for this elevation (i = 11) in Figure 6.16. Comparing the radial 
positions for these profiles enables one to see that the region where 
there is a temperature defect (negative temperatures) does not 
correspond exactly to the region where there is a flow reversal 
(negative velocities) [56,57]. The point where the velocity curve 
crosses over from positive to negative is located farther away from 
the wall. 

The physical mechanism leading to this temperature defect and 
Flow reversal will now be examined, and a short note is in order to 
‘explain how the velocity boundary layer is formed for a buoyancy- 
driven flow [58], Refer to the tvo Figures shown in Figure 6.17. In 
‘the lower figure, the total thickness of the velocity boundary layer 
can be divided into two parts. The one close to the wall is a heated 
layer where the dimensionless temperature is positive; this is the 
region where the velocity distribution is a result of balance of 
‘Suoyancy force and friction force. Moving away from the wall, there 


is an unheated layer where the dimensionless temperature is zero; the 
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Figure 6.16 Temperature and velocity profiles for NI case at i = 11 
grid Vine 
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Figure 6.17 Soundary layers on a vertical, heated surface 
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formation of the velocity in this layer is the balance of friction 
force and inertia force. 

Notice that for the buoyancy-driven flow of interest in this 
Study, the water moves ina stratified medium in which the wall 
temperature is not uniform. Then, several changes my occur when the 
water moves up toa new environnent. Refer to the top Figure in 
Figure 6.17, which is a plot of the water tenperature at a higher 
elevation. If the temperature gradient along the wall in the z 
direction is greater than that gradient of water along the 
centerline, that is dTg/dz > dT,/dz, then the thermal boundary layer 
thickness will grow. Otherwise, temperature defect will occur, 
Physically, one may say that there must be sufficient temperature 
rise in the wall in the direction of Flow in order for water to 
sustain a growing thermal boundary-layer thickness; an insufficient 
temperature rise in the wall will result in the local water 
‘temperature being lower than the ambient water at the same elevation, 
thus causing temperature defect (59-63). 

In order to determine the temperature defect region for tine 
equal to 60 seconds, the temperature gradients of the wall and water 
are plotted versus z as shown in Figure 6.18, The region where 
aT ,/dz > dT/dz is marked for temperature defect. 

The mechanism causing flow reversal can now be explained by 
again referring to Figure 6.17. In the top figure, the r position 


Where there is temperature defect can be conside 


a region of 
negative buoyancy force. This region is in the unheated layer in the 


velocity profile (lower figure) where the water velocity is mainly a 
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Figure 6.18 Determination of temperature defect region by comparing 
two temperature gradients 
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result of inertia and friction forces. Then, for water coming from 
below that is in the outer part of the velocity boundary layer, a 
negative force will arise if the inertia force in water is unable to 
overcome the negative buoyancy force above. The direction of flow 
will reverse, thus causing flow reversal. 

The fact that the crossover point of the velocity profile is 
located farther away from the wall than that of the temperature 
profile can also be explained; see Figure 6.19. Postulating that the 
entire velocity profile is the superposition of two profiles, one 
being the result of the positive branch of the temperature profile in 
Spy» and the other being the result of the temperature defect in sr, 
then the superimposed velocity profile will look like the one shown 
in solid curve in Figure 6.19(b). The crossover location for the 
velocity profile is farther away from the wall because of the 


velocity in the positive branch being larger. 


6.6 Energy Analysis 
From the thermodynamic point of view, the cooling of water in 
tank should satisfy the energy conservation principle. The initial 
energy of wall and water should be equal to the energy of the wall 
and water at any moment corrected for that portion of heat that is 
lost from the wall to the ambient air, This latter neat correction 
should be integrated over time up to that moment for energy-loss 
analysis. 
The energy analysis mentioned here serves two purposes. It 


Provides a check of the accuracy of the numerical results (other 
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Figure 6.19 Superposition of two sections of temperature profile in 
(a) to obtain velocity profiles in (>) 
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checks will be given in the next chapter). It also enables a 
comparison of the capability of the tank in preserving the energy 
stored in water. 

A plot of the difference between the initial total system energy 
and the total system energy at time t (with loss adjusted), 
normalized by the initial energy, is given in Figure 6.20. in this 
plot, a total satisfaction of the conservation principle by numerical 
results would require that this difference be zero at all time. For 
the three cases studied here, this energy conservation error ranges 
from -2 to 7% at three hours. Considering that it takes 5,400 full 
time-steps to reach three hours, and within each time-step many 
iterations are needed for computing convergent temperature, 
vorticity, and stream function values, plus the fact that there are 
approximations made in treating heat-flux continuities at interfaces, 
‘the numerical errors presented here appear reasonable. 

A comparison of the ability of the tank in preserving energy is 
made in Figure 6.21. Here the ordinate is a ratio of the energy in 
water to its initial energy, and time is plotted along the x axis. 
For all three cases studied, there is an initial rise of energy in 
water due to heat transfer from the lower wall. Al] three curves 
drop off with time because of heat loss. Of the three cases tested, 
the El case is only slightly better than II case in preserving 
energy. However, in Section 6.1.2, it has been mentioned that 
Stratification in water cannot be effectively preserved for this EI 


case. 
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Figure 6.20 Numerical errors evaluated for three cases based on the 
energy conservation principle 
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Figure 6.21 Comparison of ability of walls to preserve energy in 


water 
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6.7 Evaluation of Stratification Index 
Thermal stratification is well preserved for case II. As a 
means for further demonstration of this important result, Figure 6.22 
is given for the centerline water temperatures at 3 hours for three 
cases. The II case clearly stands out as being superior. 
A stratification index is proposed as follows to evaluate the 
Performance of storage tank in preserving stratification in a storage 


medium: 


(6.1) 


Here (@T_/dz), is the mean temperature gradient in the stratified 
Storage medium initially at the centerline of the tank. The 
differential in the numerator stands for the same temperature 
gradient at time t. ased on this definition, a perfectly stratified 
tank should have a S_ value unity at all times, and a good tank is 
one that has a high S_ value, For the three cases studied here, the 
values of St vary from 0.35 for EI, 0.46 for NI, and 0.64 for II case 
at three hours. The ability of the interior insulation to preserve 


thermal stratification is now verified quantitatively, 


Figure 6.22 Centerline water temperature profiles at 3 hours for 


three designs of the wall 
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CHAPTER VII 
REMARKS ON NUMERICAL RESULTS 


Numerical solution of a physical problem is incomplete without 
test of convergence. A test of numerical results in satisfaction of 
the energy conservation principle has been given in Chapter VI, 
Other tests must be made to show that the numerical schemes are 
stable and the numerical results are convergent, and these are the 
tests that will be made here. 

Tests of stability and convergence in this chapter will be made 
by varying the grid and time-step sizes in the computation. In the 
course of tests consideration will also be given to the computer 


(CPU) time used in the computation. 


7.1 Effect of Time-Step Size on Stability and Convergence 


Although the AOI method used in the solution of the vorticity 
transport and wall and water energy equations has been considered to 
be unconditionally stable according to the Fourier or von Neumann 
analysis, instability may still arise if boundary vorticity does not 
converge (52]. Bontoux et al. [64] made an extension of the Fourier 
Stability analysis to study the effect of boundary conditions on the 
Stability of the ADI method. They proposed a criterion for time-step 


limitation 


13a 
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1 
at z 
Ae 141 + 2oyry as 


a’ 


where both at, and ax are dimensionless. Notice that this equation 
was derived based on the analysis of natural convection in a 
differentially heated solar collector, which is, strictly speaking, 
not the kind of system being analyzed heres however, because of the 
Solution methods used in their analysis being similar to the ones 
used in the present investigation, the criterion given above can, at 
least, be served as reference for selecting time-steps for tests. 
The Prandtl number for water is taken to be 4.52 (value at 


38°C), and equation (7.1) is used to compute 


Sons om 
: 

grid size in equation (7.2) yields a time-step that is equivalent to 
areal time of 36.5 seconds. The same criterion is applied to the 
wall, where, for a 41X6 grid system, the grid size, calculated based 
on the tank radius, is 0.00533. This yields a 


‘eal time-step for the 
wall to be 4.27 seconds. This latter time-step being smaller is 
taken to be a reference for selecting test time-steps. 

Several time-step sizes (1.0, 2.0, 3.9 and 4.0 seconds) are 
tested, and the results are plotted in Figure 7.1. In Figure 7.1(a), 


the number of iterations for computing convergent neat fluxes at 
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Figure 7.1 Iteration counts for convergence tests: (a) heat flux, 
(b) boundary vorticity 
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wall/water interface is plotted versus cooling time of the tank. It 
is clear that the nunber of iterations increases rapidly with the 


time-step si 


The same trend is found when computing convergent 
boundary vorticities; see Figure 7.1(b). Here the number of 
iterations for a time-step of 4 seconds starts to oscillate, and an 
increase of time-step to § seconds does not yield convergent results, 
indicating that the numerical stability limit for the time-step has 
been exceeded. 

Since it is always desirable, if possible, to use a large time- 
step in order to reduce the number of steps computed, it is not clear 
whether such a practice is economical here because the number of 
iterations is also increased when the time-step size is increased. A 
comparison of CPU time is thus made as shown in Table 7.1. For the 
Problem under investigation, a chofce of 2 seconds time-step consumes 
‘the least computer time. 

The fluid velocities are also tested for convergence at selected 
locations and the results are presented in Figure 7.2. For the four 
time-steps tested, numerical results for three time-steps of 1.0, 
2.0, and 3.0 seconds are in close agreement at all elevations and 
radial positions tested, giving indication that the numerical results 
have converged at 2 seconds. This time-step has thus been used for 


Computing 211 numerical results presented in Chapter VI. 


7.2 Effect of Grid Size on Convergence 


In the series of tests made above, the grid system analyzed in 


water was 2 41X2I-point system. For such a system, the velocities 


Table 7.1 CPU time consumed for different 
Time steps calculated up to 
60 seconds’ cooling time 


= 

time-step 

(second) CPU time 
1.0 Lh 4nin8s 
2.0 57 min 36 s 
3.0 1h 13 min 19 5 


4.0 


233 min 34s 
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Figure 7.2 Axial velocity profiles for various time-steps at 
60 seconds 
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have shown to converge for a time-step of 2 seconds; the numerical 
computation is also stable for this time-step, while the computer 
Program is economical to run. 1t remains to be shown that this grid 
system is, in itself, good enough for a convergent solution. 

The time-step used is still taken to be 2 seconds, and one 
additional grid system (61X21) is tested by analyzing their flow and 
temperature fields in water, and the results are plotted in Figure 
7.3. A comparison of the streamlines for two grid systems shows that 
2 41X21 grid fs already able to reveal much of the flow details. On 
the other hand, the z component velocities for the 41X21 system have 
also converged; see Figure 7.4. A change of grid system to 61X21 
points in water, however, increases CPU time by a factor of 1.8 for 
60 seconds’ cooling time; an increase of grid points to 61X21 is thus 
unnecessary for the present investigation. 

A test with 41X31 grid system was also attempted, and it takes 
about 40 hours CPU time to simulate 60 seconds’ cooling time of the 
tank because of the excessive iterations needed. A prolonged 


Simulation using this grid system was not made. 


7.3. Test of Convergence Criteria for Yorticities 


It is a1s0 desirable to check if the convergence criterion for 


boundary velocities in equation ( 


38) is smal] enough to give 


accurate results. The z-component velocities are thus computed for a 
41K21 grid system with two-second time-step and three convergence 


criteria of 10°, 1074, and 10°, As shown in Figure 7.5, the 
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Figure 7.4 Axial velocity profiles for different grid systens at 
60 seconds 
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Figure 7.5 Axial velocity profiles for three convergence criteria of 
the boundary vorticity 
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velocities have converged for the 1079 criterion used in the present 


investigation. 


CHAPTER VIII 
CONCLUSIONS AND RECOMMENDATIONS 


A numerical study is presented for investigating unsteady 
natural convection in a stratified storage tank. Two different 
designs of tank wall are analyzed, one placing insulation on the 
outside of the wall, the other placing the insulation over the inside 
of the wall. An uninsulated tank is also studied for use as 
Control. For all the cases studied, water is stratified at the 
midsection of the tank; hot water occupies nearly one third at the 
top section of the tank, while cold water is located at the bottom 
section, The wall temperature is initially uniform at the hot-water 
temperature level. 

In the numerical solution, the wall and the water are treated as 
‘two Separate regions, and conditions are derived for matching 
‘temperatures and heat fluxes at the wall/water interface. External 
heat convection from the wall to the ambient air is not treated as a 
separate region, but accounted for in the boundary conditions of the 
wall, The entire problem is solved as a conjugate heat transfer 
Problem to determine the flow and temperature fields in water and the 
‘temperature field in wall. 

In the solution, the governing equations are transformed into a 
vorticity-stream function formulation, Control-volume approach is 


used to derive finite-difference equations. Kinematically consistent 
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velocity fields are used to compute vorticities and stream functions, 
and a nonuniform grid system is employed to analyze the flow and 
temperature fields in water. 

The wall and water energy equations and vorticity transport 
equation are solved by use of the ADI method, while the stream 
function equation is solved by use of the SOR method. The wall and 
water temperatures are computed simultaneously using the energy 
equations, while heat fluxes are checked at the wall/water interface 
to assure their continuity, The vorticity transport equation is 
Solved together with the stream function equation, and the 
convergence of vorticity is checked by evaluating boundary 
vorticities. Two sequences of interations are used for each full- 
time step, and a total of 5,400 full-time steps are calculated for a 
real time of three hours’ cooling of the tank. 

The performance of each storage tank is evaluated on the basis 
of its isotherm, velocity, and stream function data. Test results 
indicate that both exterior and interior insulation are about equal 
in preserving energy; however, the interior insulation performs far 
better in retaining stratification in water. Based on the 
Stratification index proposed in this dissertation, the interior 
insulation has an index value that is 85% higher than that of the 
exterior insulation. The merit of interior insulation is thus 
quantitatively established. 

The velocity and temperature in the boundary layer at the early 
Stage of cooling are also used to study flow reversal and temperature 


defect. The physical mechanisms leading to these phenomena are 


given. 


a7 


The study is concluded with a series of tests in which the 


‘time step used in the computation is shown to be numerically 


stable. 


The grid size chosen for analysis and the convergence 


criterion selected for vorticity computation are both adequate in 


yielding convergent velocities in the boundary layer. 


Based on the conclusions made above, some follow-up studies may 


be proposed as follows: 


a 


(2) 


(3) 


The study made here is based on one aspect ratio. It is 
desirable to test other aspect ratios so that an optimum 
configuration of the tank may be determined. 

in the present study, the performance of the tank is 
analyzed under a static condition. It may be necessary to 
evaluate the performance of the tank under a dynamic 
condition. 

The numerical study presented here can be made more useful 
if it is verified by experiment. A limited experimental 
Program serves well to meet this need. Once verified to be 
accurate, the computer program developed here can be used to 
test other operating conditions, thereby leading to an 


optimum use of the storage tank. 


APPENDIX A 


PHYSICAL PROPERTIES: 


= 
water air 1 carbon | soft 
(38°C) | (10°C) steel rubber 
specific heat 4.187 0.473 1.884 
ka/Kg"C 
density 993.1 7800 10 
by kg/m 
conductivity 0.83 0.028 43.0 0.173 
ke win" | 
L 
coefficient of a.exiom | 3,53x10°S | 
expansion | 
| 
8, Ko | 
kinematic viscosity | 6.87K1077 | 1, 36x10 
vy m/s 
diffusivity 1,576x1077 aetoxio"S | 827x107) 
a, m/s 
Prandtl nunber 4.52 ona 
Pr 
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APPENDIX B 
TIME SCALE ANALYSIS FOR DETERMINING 
HEAT TRANSFER MECHANISMS 


The time-scale analysis given in this Appendix provides a 
Qualitative prediction of the preferred heat transfer mechanisms in a 
storage tank. Such an analysis has also been made by others 
(58,65,66] to study other natural convection problems. 


Consider the energy equation in the tank wall 


(3.24) 
and the energy equation in the water medium 
ar a, aT, fie pate oF, 
fev ate, he atte tot, of (az) 
Fe a ae toe i 


Six time scales can be derived from these tvo equations and these 
time scales can be used to relate six heat transfer mechanisms as 


‘illustrated in Figure 8.1. These time scales are defined as follows: 
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wall water 
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I 
Sf.b | 
| 
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| boundary stratified 
Tayer core 
/ resion region 
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Figure 8.1 Six time scales for six heat-transfer mechanisms 


Asi 


ts,2? time scale for z-direction heat conduction in the wall, 
typi time scale for r-direction heat conduction in the wall, 


te,q time scale for r-direction heat conduction across the 
“boundary layer, 


tf,pi tine scale for z-direction heat convection in the boundary 
0" layer, 


tr,2! time scale for z-direction heat diffusion in the core 
"7 region, 


tp,ci time scale for z-direction heat convection in the core 
°°" region, 
‘An order-of-nagni tude analysis can be used to derive the 


expressions for time scales as follows. 


Consider the energy equation (3.24) for heat conduction in the z 
direction, Uropping the radial differential terms in this equation 


gives 


fis, (3.1) 
¢ 


where FH is the fraction of the height of the wall along which there 
is a large temperature gradient; H is the height of this wall. This 


equation can be rearranged to be 
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t, + (3.2) 


A Schematic diagram showing the wall length over which the 


temperature change is al is provided in Figure 8.2(a). 
(i) ty pe 


Consider equation (3.24) in the r direction. Dropping the z 


differential term in this equation gives 


wal gah (8.3) 
Bp oe 


Since Rj>>8, the second term in the parenthesis is small. It follows 


that 


2 
fo ai (3.4) 


Sr 
s 


The temperature gradient in the radial direction in the wall is shown 


schematically in Figure 8.2(b). 
GH teat 
For heat transfer across the boundary layer, sy, the temperature 


drop across this layer is aT; see Figure 8.2(c). Oropping the 


convection terms and the z-direction conduction term in equation 


abe 


(a) (b) 


Fat 


in 
aie 
— 


Figure 8,2 Schematic diagrams for studying time scales 
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(3.12) gives 


(3.5) 
This equation can be simplified to be 
2 
& 
te, tae (3.6) 
fan "ag 


Gv) tpt 


For heat convection in the boundary layer, those terms on the 


right-hand side of equation (3.12) become negligibly small, giving 


tL 
Fb 


at 
~ Yau For (8.7) 


where Fal represents the fraction of the flow length along the wall 
where the buoyancy force has an effect on the Fluid motion; see 
Figure 8.2(d), The quantity Vz 9, is used to represent the 2 
velocity inside the boundary layer. Equation (8.7) can be arranged 


to write 


te te (3.8) 


iss 


(Wy) thet 


For the diffusion of heat in the core region, the convection 


terms in equation (3.12) vanish, giving 


zy, 


fu 
a, (3.9) 

fpf UF a) 

where Fell represents the fraction of the tank height along which the 

diffusion of heat is significant in the core region; see Figure 

3.2(e). A rearrangement of equation (8.9) leads to 


242 
Few 


te * (8,10) 


For convective flow in the core region, the diffusion terms on 


the right-hand side of equation (3.12) vanish, yielding 


fi or 
“¥, (a.1) 
The Ft 


where Vz represents the z velocity in the core region; see Figure 
3.2(f). It 1s anticipated that, in the core region, the displacement 
of water in the z direction is closely related to the convective 
sotion of water in the boundary region; Fi can thus be used to 


Tepresent the displacement in the core region. The time scale for 
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convection in the core region can be derived by rearranging equation 


(8.11) as 


FH 
the (3,12) 


A close examination of the six time-scale expressions derived 
above reveals that there are three unknowns (S;, Vp gu» Vz) in these 
expressions. These three unknowns oust be expressed in terms of 
known quantities before the time scales can be calculated. Use is 
thus made of the continuity equation, the momentum equation, and the 


energy equation in the boundary layer 


(8.13) 


Z + 98(T - T,) (3.14) 


(3,15) 


In these equations, the length scales r and z may be expressed as 6 
and Fy, respectively, As expected, Fy is much greater than &. 
The continuity equation (8,13) can be used to derive 


B 


ae (8.16) 
CaS 


From the energy equation (3.15), one can derive two relations 


git 
ree 


From equation (8.17) gives 


ar H 
vod 


T 
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(a.a7) 


(8.18) 


(a.19) 


The momentum equation (8.14) can be used to balance the viscous 


that 


v ay FH 

Pee as 3 oy 
geet ~ vy Eg = 
T T 


A slight rearrangement of this equation gives 


t 1 
avr é 
faa 7 
+ (Rp) + Fal a 

3 


force and the bouyancy force in the boundary layer [58]. 


It follows 


(8.20) 


(8.21) 


From conservation of mss between the boundary layer flow and 


‘the flow in the core region gives 


iz 


6 
i 
MeO Ry Yaya 


(3,22) 
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The time-scale expressions given above are now used to study a 
sample case in which the storage tank is not insulated; it is desired 
to compare the importance of various heat transfer mechanisms at 60 
seconds. The isotherms plotted for this tine (Figure 6.1(a)) are 
used to evaluate 


1 
eed 
(8.23) 
1 
ees 
Other data are listed as follows: 
“5 42 
ty 7 1.1610 a ds 
"* 1576x1077 ws 
a.27x08 a/s 
(3.24) 


6.87x107 a/s 
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With the use of these data for input in equations (8.19), (8.21) 


and (8.22), values for ér, Vp a1, and Vp can be estimated as 


4 = 8.a5K0™ @ 
-2 

¥, au + 4:02x10 n/s (3.25) 

v, = 2.3710 a/s 


In contrast, the computed results are 


4, = 58 of Ry = 7.5x10 


. 


2a, = 3+8X10" m/s (maximum) (8.27) 


5.6x10"* m/s (maximum) 


There is a difference of one order of magnitude between the estimated 

& and the computed 67, and also between the estinated Vz,q) and the 

computed Vz 91+ However, the estimation of ¥, appears satisfactory. 
Introducing the estimated values of 57, Vp gi» and Yq into the 


time-scale expressions gives 


ty,2 7 7-76K10" seconds t, 2 = 1.09x10° seconds 
typ = 1638 seconds tp = 193.47 seconds 
tpg = 4:97 seconds 


(8,28) 
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ty) = 4497 seconds 


fs 
t, =2, 54x0° seconds 
faa 


2 
Teg = 844x107 seconds 


Some comments can be made based on the time-scale values calculated 


here, 


CO) ty petyys 


This relation implies that conduction in the z-direction in the 
wall is much slower than conduction in the r direction in the wall. 
‘The temperature distribution in the wall is primarily determined by 


‘the conduction in the z direction. 
tha * te 


This relation implies that conduction and convection are equally 


important in the boundary layer. 


(iH) ty cete,» 


This relation implies that in the core region heat diffusion is 
‘much Slower than heat convection; only convection is important in the 


core region. 
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(iv) ty,28te,2 and ty,2/te,¢ = 9.19 


Heat diffusion in the core region is a much slower process than 
heat conduction in the wall in the z direction. Heat convection in 
the core region, however, is only about an order of magnitude faster 


‘than the heat conduction in the wall in the z direction. 


) tealtese * teyp/te cl 


Heat convection in the core region is a much slower process than 
either the heat conduction across, or the heat convection in, the 
boundary Tayer. Heat transfer in the boundary layer is virtually 
unaffected by tne heat transfer in the core region. 

The analysis made above can also be extended to study heat 


‘transfer in an insulated tank. 


(VI) ty pst neta 


Heat transfer across the thickness of the insulation is much 
Slower than either conduction across the tank wall or conduction 
across the boundary layer. Heat transfer in the radial direction is 


‘thus controlled mainly by the thermal resistance in the insulation. 
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(Wi) tyy2tety, 2 and ty,7 # lyn + 2b pk rete 


In the z direction heat conduction in the insulation is much 
Slower than heat conduction in the wall. The second relation 
signifies that, if the insulation is placed inside, it is much faster 
for the heat transfer to take the route of conduction in the z 
direction fn the wall than either conduction in the z direction in 


‘the insulation or diffusion of heat in the core region. 


WAH) typ * th gta 


Here the relation implies that, when the insulation is placed 
outside, the time required for the heat to transfer across the 
boundary layer and wall is much faster than the heat to transfer 
across the insulation. Using this relation together with the first 
relation given in (vii) above, one may conclude that the heat will 
leak more conveniently through the tank wall than through the 


insulation outside. 
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